home *** CD-ROM | disk | FTP | other *** search
/ Aminet 6 / Aminet 6 - June 1995.iso / Aminet / gfx / 3d / irit50src.lha / irit5 / symb_lib / arc_len.c < prev    next >
Encoding:
C/C++ Source or Header  |  1994-11-09  |  17.0 KB  |  443 lines
  1. /******************************************************************************
  2. * Arc_len.c - arc length functions and unit length normalization.          *
  3. *******************************************************************************
  4. * Written by Gershon Elber, Apr. 93.                          *
  5. ******************************************************************************/
  6.  
  7. #include "symb_loc.h"
  8.  
  9. #define MAX_ALEN_IMPROVE_ITERS 5
  10. #define ARCLEN_CONST_SET_EPSILON 0.001
  11.  
  12. /*****************************************************************************
  13. * DESCRIPTION:                                                               M
  14. * Subdivides the given curves to curves, each with size of control polygon   M
  15. * less than or equal to MaxLen.    Returned is a list of curves.             M
  16. *                                                                            *
  17. * PARAMETERS:                                                                M
  18. *   Crv:        To subdivide into curves, each with control polygon length   M
  19. *               less than MaxLen.                         M
  20. *   MaxLen:     Maximum length of control polygon to allow.                  M
  21. *                                                                            *
  22. * RETURN VALUE:                                                              M
  23. *   CagdCrvStruct *:   List of subdivided curves from Crv, each with control M
  24. *                      polygon size of less than MaxLen.                     M
  25. *                                                                            *
  26. * KEYWORDS:                                                                  M
  27. *   SymbLimitCrvArcLen, arc length                                           M
  28. *****************************************************************************/
  29. CagdCrvStruct *SymbLimitCrvArcLen(CagdCrvStruct *Crv, CagdRType MaxLen)
  30. {
  31.     if (SymbCrvArcLenPoly(Crv) > MaxLen) {
  32.     CagdRType TMin, TMax;
  33.     CagdCrvStruct *Crv1, *Crv2, *Crv1MaxLen, *Crv2MaxLen;
  34.  
  35.     CagdCrvDomain(Crv, &TMin, &TMax);
  36.  
  37.     Crv1 = CagdCrvSubdivAtParam(Crv, (TMin + TMax) / 2.0);
  38.     Crv2 = Crv1 -> Pnext;
  39.  
  40.     Crv1MaxLen = SymbLimitCrvArcLen(Crv1, MaxLen);
  41.     Crv2MaxLen = SymbLimitCrvArcLen(Crv2, MaxLen);
  42.  
  43.     CagdCrvFree(Crv1);
  44.     CagdCrvFree(Crv2);
  45.  
  46.     for (Crv1 = Crv1MaxLen; Crv1 -> Pnext != NULL; Crv1 = Crv1 -> Pnext);
  47.     Crv1 -> Pnext = Crv2MaxLen;
  48.     return Crv1MaxLen;
  49.     }
  50.     else
  51.         return CagdCrvCopy(Crv);
  52. }
  53.  
  54. /*****************************************************************************
  55. * DESCRIPTION:                                                               M
  56. * Computes a bound on the arc length of a curve by computing the length of   M
  57. * its control polygon.                                 M
  58. *                                                                            *
  59. * PARAMETERS:                                                                M
  60. *   Crv:       To bound its length.                                          M
  61. *                                                                            *
  62. * RETURN VALUE:                                                              M
  63. *   CagdRType:   An upper bound on the curve Crv length as the length of     M
  64. *                Crv's control polygon.                                      M
  65. *                                                                            *
  66. * KEYWORDS:                                                                  M
  67. *   SymbCrvArcLenPoly, arc length                                            M
  68. *****************************************************************************/
  69. CagdRType SymbCrvArcLenPoly(CagdCrvStruct *Crv)
  70. {
  71.     int i;
  72.     CagdCrvStruct
  73.     *CrvE3 = CagdCoerceCrvTo(Crv, CAGD_PT_E3_TYPE);
  74.     CagdRType
  75.     Len = 0.0,
  76.     **Points = CrvE3 -> Points;    
  77.  
  78.     for (i = 1; i < CrvE3 -> Length; i++)
  79.         Len += sqrt(SQR(Points[1][i] - Points[1][i - 1]) +
  80.             SQR(Points[2][i] - Points[2][i - 1]) +
  81.             SQR(Points[3][i] - Points[3][i - 1]));
  82.  
  83.     CagdCrvFree(CrvE3);
  84.  
  85.     return Len;
  86. }
  87.  
  88. /*****************************************************************************
  89. * DESCRIPTION:                                                               M
  90. * Normalizes the given vector field curve to be a unit length curve, by      M
  91. * computing a scalar curve to multiply with this vector field curve.         M
  92. *   Returns the multiplied curve if Mult, or otherwise just the scalar       M
  93. * curve.                                     M
  94. *                                                                            *
  95. * PARAMETERS:                                                                M
  96. *   OrigCrv:     Curve to approximate a unit size for.                 M
  97. *   Mult:        Do we want to multiply the computed scalar curve with Crv?  M
  98. *   Epsilon:     Accuracy required of this approximation.                    M
  99. *                                                                            *
  100. * RETURN VALUE:                                                              M
  101. *   CagdCrvStruct *:  A scalar curve to multiply OrigCrv so a unit size      M
  102. *                     curve will return if Mult is FALSE, or the actual      M
  103. *                     unit size vector field curve, if Mult.             M
  104. *                                                                            *
  105. * KEYWORDS:                                                                  M
  106. *   SymbCrvUnitLenScalar, unit vector field                                  M
  107. *****************************************************************************/
  108. CagdCrvStruct *SymbCrvUnitLenScalar(CagdCrvStruct *OrigCrv,
  109.                     CagdBType Mult,
  110.                     CagdRType Epsilon)
  111. {
  112.     int i, j;
  113.     CagdCrvStruct
  114.     *ScalarCrv = NULL,
  115.     *Crv = CAGD_IS_BEZIER_CRV(OrigCrv) ? CnvrtBezier2BsplineCrv(OrigCrv)
  116.                        : CagdCrvCopy(OrigCrv);
  117.     CagdBType
  118.         IsRational = CAGD_IS_RATIONAL_CRV(Crv);
  119.  
  120.     for (i = 0; i < MAX_ALEN_IMPROVE_ITERS; i++) {
  121.     CagdCrvStruct *UnitCrvAprx, *ScalarCrvSqr,
  122.         *DotProdCrv = SymbCrvDotProd(Crv, Crv);
  123.     CagdRType Min, Max, *SCPoints,
  124.         *DPPoints = DotProdCrv -> Points[1];
  125.  
  126.     if (ScalarCrv)
  127.         CagdCrvFree(ScalarCrv);
  128.     ScalarCrv = CagdCrvCopy(DotProdCrv);
  129.     SCPoints = ScalarCrv -> Points[1];
  130.  
  131.     /* Compute an approximation to the inverse function. */
  132.     for (j = 0; j < ScalarCrv -> Length; j++, DPPoints++) {
  133.         *SCPoints++ = *DPPoints > 0.0 ? 1.0 / sqrt(*DPPoints) : 1.0;
  134.     }
  135.     ScalarCrvSqr = SymbCrvMult(ScalarCrv, ScalarCrv);
  136.  
  137.     /* Multiply the two to test the error. */
  138.     UnitCrvAprx = SymbCrvMult(ScalarCrvSqr, DotProdCrv);
  139.     CagdCrvFree(ScalarCrvSqr);
  140.  
  141.     CagdCrvMinMax(UnitCrvAprx, 1, &Min, &Max);
  142.  
  143.     if (1.0 - Min < Epsilon && Max - 1.0 < Epsilon) {
  144.         CagdCrvFree(UnitCrvAprx);
  145.         CagdCrvFree(DotProdCrv);
  146.         break;
  147.     }
  148.     else {
  149.         /* Refine in regions where the error is too large. */
  150.         int k,
  151.             Length = UnitCrvAprx -> Length,
  152.             Order = UnitCrvAprx -> Order,
  153.             KVLen = Length + Order;
  154.         CagdRType
  155.         *KV = UnitCrvAprx -> KnotVector,
  156.         *RefKV = IritMalloc(sizeof(CagdRType) * 2 * Length),
  157.             *Nodes = BspKnotNodes(KV, KVLen, Order),
  158.         **Points = UnitCrvAprx -> Points;
  159.  
  160.         for (j = k = 0; j < Length; j++) {
  161.         CagdRType
  162.             V = 1.0 - (IsRational ? Points[1][j] / Points[0][j]
  163.                       : Points[1][j]);
  164.  
  165.         if (ABS(V) > Epsilon) {
  166.             int Index = BspKnotLastIndexLE(KV, KVLen, Nodes[j]);
  167.  
  168.             if (APX_EQ(KV[Index], Nodes[j])) {
  169.             if (j > 0)
  170.                 RefKV[k++] = (Nodes[j] + Nodes[j - 1]) / 2.0;
  171.             if (j < Length - 1)
  172.                 RefKV[k++] = (Nodes[j] + Nodes[j + 1]) / 2.0;
  173.             }
  174.             else
  175.             RefKV[k++] = Nodes[j];
  176.         }
  177.         }
  178.         CagdCrvFree(UnitCrvAprx);
  179.         CagdCrvFree(DotProdCrv);
  180.  
  181.         IritFree((VoidPtr) Nodes);
  182.  
  183.         if (k == 0) {
  184.         /* No refinement was found needed - return current curve. */
  185.         IritFree((VoidPtr) RefKV);
  186.         break;
  187.         }
  188.         else {
  189.         CagdCrvStruct
  190.             *CTmp = CagdCrvRefineAtParams(Crv, FALSE, RefKV, k);
  191.  
  192.         IritFree((VoidPtr) RefKV);
  193.         CagdCrvFree(Crv);
  194.         Crv = CTmp;
  195.         }
  196.  
  197.     }
  198.     }
  199.  
  200.     CagdCrvFree(Crv);
  201.  
  202.     /* Error is probably within bounds - returns this unit length approx. */
  203.     if (Mult) {
  204.     CagdCrvStruct *CrvX, *CrvY, *CrvZ, *CrvW, *CTmp;
  205.     int MaxCoord = CAGD_NUM_OF_PT_COORD(OrigCrv -> PType);
  206.  
  207.     SymbCrvSplitScalar(ScalarCrv, &CrvW, &CrvX, &CrvY, &CrvZ);
  208.     CagdCrvFree(ScalarCrv);
  209.     ScalarCrv = SymbCrvMergeScalar(CrvW,
  210.                        CrvX,
  211.                        MaxCoord > 1 ? CrvX : NULL,
  212.                        MaxCoord > 2 ? CrvX : NULL);
  213.     CagdCrvFree(CrvX);
  214.     if (CrvW)
  215.         CagdCrvFree(CrvW);
  216.  
  217.     CTmp = SymbCrvMult(ScalarCrv, OrigCrv);
  218.     CagdCrvFree(ScalarCrv);
  219.  
  220.     return CTmp;
  221.     }
  222.     else {
  223.     return ScalarCrv;
  224.     }
  225. }
  226.  
  227. /*****************************************************************************
  228. * DESCRIPTION:                                                               M
  229. * Computes the curve which is a square root approximation to a given scalar  M
  230. * curve, to within epsilon.                             M
  231. *                                                                            *
  232. * PARAMETERS:                                                                M
  233. *   OrigCrv:    Scalar curve to approximate its square root function.        M
  234. *   Epsilon:    Accuracy of approximation.                                   M
  235. *                                                                            *
  236. * RETURN VALUE:                                                              M
  237. *   CagdCrvStruct *:   A curve approximating the square root of OrigCrv.     M
  238. *                                                                            *
  239. * KEYWORDS:                                                                  M
  240. *   SymbCrvSqrtScalar, square root                                           M
  241. *****************************************************************************/
  242. CagdCrvStruct *SymbCrvSqrtScalar(CagdCrvStruct *OrigCrv, CagdRType Epsilon)
  243. {
  244.     int i, j;
  245.     CagdCrvStruct
  246.     *ScalarCrv = NULL,
  247.     *Crv = CAGD_IS_BEZIER_CRV(OrigCrv) ? CnvrtBezier2BsplineCrv(OrigCrv)
  248.                        : CagdCrvCopy(OrigCrv);
  249.     CagdBType
  250.         IsRational = CAGD_IS_RATIONAL_CRV(Crv);
  251.  
  252.     for (i = 0; i < MAX_ALEN_IMPROVE_ITERS; i++) {
  253.     CagdCrvStruct *ScalarCrvSqr, *ErrorCrv;
  254.     CagdRType Min, Max, *SCPoints,
  255.         *Points = Crv -> Points[1],
  256.         *WPoints = IsRational ? Crv -> Points[0] : NULL;
  257.  
  258.     if (ScalarCrv)
  259.         CagdCrvFree(ScalarCrv);
  260.     ScalarCrv = CagdCrvCopy(Crv);
  261.     SCPoints = ScalarCrv -> Points[1];
  262.  
  263.     /* Compute an approximation to the inverse function. */
  264.     for (j = 0; j < ScalarCrv -> Length; j++) {
  265.         CagdRType
  266.         V = IsRational ? *Points++ / *WPoints++ : *Points++;
  267.  
  268.         *SCPoints++ = V >= 0.0 ? sqrt(V) : 0.0;
  269.     }
  270.     ScalarCrvSqr = SymbCrvMult(ScalarCrv, ScalarCrv);
  271.  
  272.     /* Compare the two to test the error. */
  273.     ErrorCrv = SymbCrvSub(ScalarCrvSqr, Crv);
  274.     CagdCrvFree(ScalarCrvSqr);
  275.  
  276.     CagdCrvMinMax(ErrorCrv, 1, &Min, &Max);
  277.  
  278.     if (Min > -Epsilon && Max < Epsilon) {
  279.         CagdCrvFree(ErrorCrv);
  280.         break;
  281.     }
  282.     else {
  283.         /* Refine in regions where the error is too large. */
  284.         int k,
  285.             Length = ErrorCrv -> Length,
  286.             Order = ErrorCrv -> Order,
  287.             KVLen = Length + Order;
  288.         CagdRType
  289.         *KV = ErrorCrv -> KnotVector,
  290.         *RefKV = IritMalloc(sizeof(CagdRType) * 2 * Length),
  291.             *Nodes = BspKnotNodes(KV, KVLen, Order),
  292.         **ErrPoints = ErrorCrv -> Points;
  293.  
  294.         for (j = k = 0; j < Length; j++) {
  295.         CagdRType
  296.             V = 1.0 - (IsRational ? ErrPoints[1][j] / ErrPoints[0][j]
  297.                       : ErrPoints[1][j]);
  298.  
  299.         if (ABS(V) > Epsilon) {
  300.             int Index = BspKnotLastIndexLE(KV, KVLen, Nodes[j]);
  301.  
  302.             if (APX_EQ(KV[Index], Nodes[j])) {
  303.             if (j > 0)
  304.                 RefKV[k++] = (Nodes[j] + Nodes[j - 1]) / 2.0;
  305.             if (j < Length - 1)
  306.                 RefKV[k++] = (Nodes[j] + Nodes[j + 1]) / 2.0;
  307.             }
  308.             else
  309.             RefKV[k++] = Nodes[j];
  310.         }
  311.         }
  312.         CagdCrvFree(ErrorCrv);
  313.  
  314.         IritFree((VoidPtr) Nodes);
  315.  
  316.         if (k == 0) {
  317.         /* No refinement was found needed - return current curve. */
  318.         IritFree((VoidPtr) RefKV);
  319.         break;
  320.         }
  321.         else {
  322.         CagdCrvStruct
  323.             *CTmp = CagdCrvRefineAtParams(Crv, FALSE, RefKV, k);
  324.  
  325.         IritFree((VoidPtr) RefKV);
  326.         CagdCrvFree(Crv);
  327.         Crv = CTmp;
  328.         }
  329.  
  330.     }
  331.     }
  332.  
  333.     CagdCrvFree(Crv);
  334.  
  335.     return ScalarCrv;
  336. }
  337.  
  338. /*****************************************************************************
  339. * DESCRIPTION:                                                               M
  340. * Computes a scalar curve approximating the arc length of given curve Crv.   M
  341. * Arc length is astimated by computing the square of Crv first derivative    M
  342. * approximating its square root and integrating symbolically.             M
  343. *                                                                            *
  344. * PARAMETERS:                                                                M
  345. *   Crv:       To approximate its arc length scalar field.                   M
  346. *   Epsilon:   Accuracy of approximating.                                    M
  347. *                                                                            *
  348. * RETURN VALUE:                                                              M
  349. *   CagdCrvStruct *:   A scalar field approximating Crv arc length.          M
  350. *                                                                            *
  351. * KEYWORDS:                                                                  M
  352. *   SymbCrvArcLenCrv, arc length                                             M
  353. *****************************************************************************/
  354. CagdCrvStruct *SymbCrvArcLenCrv(CagdCrvStruct *Crv, CagdRType Epsilon)
  355. {
  356.     CagdCrvStruct
  357.     *DerivCrv = CagdCrvDerive(Crv),
  358.         *DerivMagSqrCrv = SymbCrvDotProd(DerivCrv, DerivCrv),
  359.     *DerivMagCrv = SymbCrvSqrtScalar(DerivMagSqrCrv, Epsilon),
  360.         *ArcLenCrv = CagdCrvIntegrate(DerivMagCrv);
  361.  
  362.     CagdCrvFree(DerivCrv);
  363.     CagdCrvFree(DerivMagSqrCrv);
  364.     CagdCrvFree(DerivMagCrv);
  365.  
  366.     return ArcLenCrv;
  367. }
  368.  
  369. /*****************************************************************************
  370. * DESCRIPTION:                                                               M
  371. * Computes a tight approximation to the arc length of a curve.             M
  372. *   Estimates the arc length scalar field of Crv using SymbCrvArcLenCrv and  M
  373. * evaluate the estimate on the curve's domain boundary.                      M
  374. *                                                                            *
  375. * PARAMETERS:                                                                M
  376. *   Crv:        Curve to compute a tight approximation on arc length.        M
  377. *   Epsilon:    Accuracy control.                                            M
  378. *                                                                            *
  379. * RETURN VALUE:                                                              M
  380. *   CagdRType:   The approximated arc length of the given curve Crv.         M
  381. *                                                                            *
  382. * KEYWORDS:                                                                  M
  383. *   SymbCrvArcLen, arc length                                                M
  384. *****************************************************************************/
  385. CagdRType SymbCrvArcLen(CagdCrvStruct *Crv, CagdRType Epsilon)
  386. {
  387.     CagdRType TMin, TMax, *Pt;
  388.     CagdCrvStruct
  389.     *ArcLenCrv = SymbCrvArcLenCrv(Crv, Epsilon);
  390.  
  391.     CagdCrvDomain(ArcLenCrv, &TMin, &TMax);
  392.     Pt = CagdCrvEval(ArcLenCrv, TMax);
  393.     CagdCrvFree(ArcLenCrv);
  394.  
  395.     return CAGD_IS_RATIONAL_CRV(ArcLenCrv) ? Pt[1] / Pt[0] : Pt[1];
  396. }
  397.  
  398. /*****************************************************************************
  399. * DESCRIPTION:                                                               M
  400. * Computes parameter values to move steps of Length at a time on curve Crv.  M
  401. *    Returned is a list of parameter values to move along.             M
  402. *                                                                            *
  403. * PARAMETERS:                                                                M
  404. *   Crv:        Curve to compute constant arc Length steps.                  M
  405. *   Length:     The step size.                                               M
  406. *   Epsilon:    Accuracy control.                                            M
  407. *                                                                            *
  408. * RETURN VALUE:                                                              M
  409. *   CagdPtStruct *:  List of parameter values to march along Crv with arc    M
  410. *                    Length between them.                                    M
  411. *                                                                            *
  412. * KEYWORDS:                                                                  M
  413. *   SymbCrvArcLenSteps, arc length                                           M
  414. *****************************************************************************/
  415. CagdPtStruct *SymbCrvArcLenSteps(CagdCrvStruct *Crv,
  416.                  CagdRType Length,
  417.                  CagdRType Epsilon)
  418. {
  419.     CagdRType TMin, TMax, *Pt, CrvLength, Len;
  420.     CagdPtStruct *Param,
  421.     *ParamList = NULL;
  422.     CagdCrvStruct
  423.     *ArcLenCrv = SymbCrvArcLenCrv(Crv, Epsilon);
  424.  
  425.     CagdCrvDomain(ArcLenCrv, &TMin, &TMax);
  426.     Pt = CagdCrvEval(ArcLenCrv, TMax);
  427.  
  428.     CrvLength = CAGD_IS_RATIONAL_CRV(ArcLenCrv) ? Pt[1] / Pt[0] : Pt[1];
  429.  
  430.     for (Len = CrvLength - Length; Len > 0.0; Len -= Length) {
  431.     Param = SymbCrvConstSet(ArcLenCrv, 1, ARCLEN_CONST_SET_EPSILON, Len);
  432.     if (Param == NULL || Param -> Pnext != NULL)
  433.         SYMB_FATAL_ERROR(SYMB_ERR_REPARAM_NOT_MONOTONE);
  434.  
  435.     Param -> Pnext = ParamList;
  436.     ParamList = Param;
  437.     }
  438.  
  439.     CagdCrvFree(ArcLenCrv);
  440.  
  441.     return ParamList;
  442. }
  443.